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(54) Abstract Title 

Filtering noise from seismic data 



(57) A method relating to filtering coherent noise and interference from seismic data by constrained adaptive 
beamforming is described using a constraint design methodology which allows the imposition of an arbitrary 
predesigned quiescent response on the beamformer. The method also makes sure that the beamformer 
response in selected regions of the frequency-wavenumber space is entirely controlled by this quiescent 
response, hence ensuring signal preservation and robustness to perturbations. Built-in regularization brings 
an additional degree of robustness. Seismic signals with arbitrary spectral content in the 
frequency-wavenumber domain are preserved, while coherent noise and interference that is temporally and 
spatially nonstationary is adaptively filtered. The approach is applicable to attenuation of all types of coherent 
noise in seismic data including swell-noise, bulge-wave noise, ground-roll, airwave, seismic vessel and rig 
interference, etc. It is applicable to both linear or areal arrays. 
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specification of the desired quiescent weight vector, 
which defines the first column of the constraint matrix 



specification of the signal protection region A 
( in the frequency-wavenumber space) 



computation of the correlation matrix 
of all the projected steering vectors in region A 



determination of the principal eigenvectors of thecorrelation 
matrix as the remaining columns of the constraint matrix 



specification of the of the constraint matrix 

and the response vector 



FIG. 4 
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Performing seismic acquisition to acquire single sensor 
recordings representing seismic waves energy 



filtering single sensor recordings prior to group forming 
using a ( partially ) adaptive filter with constraints 



generating single sensor recordings 
with reduced amount of noise 



further processing the single sensor recordings 
with reduced amount of noise 



generating a representation of subterranean formation(s) 



FIG. 4 
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" 1 " 

Adaptive Seismic Noise and Interference Attenuation Method 

This invention relates to seismic data acquisition and to 
methods of processing seismic data. It relates to a process for 
filtering coherent noise and interference from seismic data 
by an adaptive beamforming method. In another aspect, it relates 
5 to adaptively filtering coherent noise and interference from 
seismic data while preserving seismic signals with arbitrary 
spectral content in the frequency- wavenumber domain. In yet 
another aspect, it relates to adaptively filtering coherent 
noise and interference that is temporally and spatially 
10 nonstationary . In a further aspect, it relates to adaptively 

filtering coherent noise and interference that has been recorded 
by a sensor array in the presence of perturbations. 

15 BACKGROUND OF THE INVENTION 

In seismic surveys, a seismic source induces seismic waves at or 
near the surface of the earth. Explosives, vibrating devices and 
airguns are examples of seismic sources. These waves propagate 

20 through the earth and are reflected, refracted, and diffracted 
by formations within the earth, and can be detected by a 
plurality of sensors (typically geophones or hydrophones) at the 
earth's surface. Each such receiver monitors the seismic 
wavefield, which is then recorded. The data received by a 

25 receiver and then recorded are collectively called a trace. The 
collection of traces are stored for further processing to gain 
information about the earth's subsurface. Such information is 
commonly interpreted to detect the possible presence of 
hydrocarbons, or to monitor changes in hydrocarbon bearing 

30 rocks. . . 
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Seismic data in general contains coherent noise signals, along 
with seismic reflection signals. These noise signals, hereafter 
referred to as the noise, interfere with the interpretation of 
the seismic signals, and degrade the quality of the subsurface 
5 images that can be obtained by further processing. It is 

therefore desirable to suppress the noise that is present in 
the recorded data before processing it for imaging. 

In land seisrnics, source-generated noise like ground- roll and 
10 air-waves are the dominant noise types, and can lead to severe 
degradation in data quality. In marine seisrnics, energy 
propagating as waves trapped in the water column and near- 
surface layers is a significant source, as well as swell noise 
and bulge -wave noise which result from waves propagating 
15 along the streamers of receiver devices. Other sources of 

coherent noise in marine seisrnics include passing vessels, other 
vessels acquiring seismic data in the vicinity, or any 
drilling activity close to the survey area. 

20 An important feature of the so-called coherent noise present in 
seismic data is the distance over which the noise appears 
coherent. In many circumstances, the noise is coherent over 
only a few meters. In other cases, although the noise is mostly 
coherent, there exists spatially impulse noise. In such cases, 

25 filtering methods which have large spatial extent, like the 
known f requency-wavenumber filtering generate undesirable 
artifacts, which are mistakenly identified as seismic events 
after further processing and imaging. 

30 Another feature of the noise present in seismic data is that it 
is often non- stationary as a function of time; i.e. its 
characteristics change as a function of time. 
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During recent years there have been suggested a variety of 
methods employing the central concept of applying adaptive 
signal processing ideas to the problem of suppressing coherent 
noise in seismic data. 

5 

Booker and Ong, in: "Multiple-constraint adaptive filtering," 
Geophysics, Vol.36, pp. 498-509, 1971, have derived an algorithm 
for multi-channel time-series data processing, which maintains 
specified initial multiple filter constraints for known signal 
10 or noise sources while simultaneously adapting the filter to 

minimize the effect of the unknown source field. The constraints 
are of the multiple look-direction constraints type, where look- 
directions must be precisely specified. 



15 In the International Patent Application W097/25632, Ozbek has 
described a class of adaptive signal processing techniques for 
attenuation of dispersive, nonstationary and aliased coherent 
noise in seismic data, in the presence of phase and amplitude 
perturbations. The methods developed can be classified as multi- 

20 channel adaptive interference cancellers. Since a signal -free 
noise reference is not readily available in seismic data 
acquisition, various preprocessing techniques are introduced to 
generate the coherent noise reference channels. In the single- 
component version of the method, moveout (apparent velocity) and 

25 spatio-temporal coherence are used as the criteria for 

differentiating between the signal and the noise. In the multi- 
component version, polarization is used as an additional 
attribute for differentiation. Once single or multiple 
noise reference channels are established, coherent noise in the 

30 primary channel is canceled using data-adaptive least-squares 
multi -channel filter banks. 
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US Patent No. 4,556,962 attempts to attenuate the ground roll 
from a surface seismic source by placing a sensor close to the 
source to detect the interfering noise. The interfering noise 
is scaled, delayed and summed with signals from a more distant 
5 geophone array and then cross -correlated with the original 

vibrational source. This patent also suggests that an adaptive 
filter be used so as to modify the delayed signal to correspond 
more closely to that detected by the more distant geophone 
array. However, ground roll is in general of an inhomogeneous 

10 nature; due to dispersion and scattering from near surface 
anomalies the ground roll measured at one point increasingly 
deviates in character from that measured at another with 
increasing distance. Hence, the ground roll measured close to 
the source may be substantially different from that received by 

15 the geophone array, and the adaptive filter may not be able 
to deal with this. It is also difficult to measure seismic 
signals (ground roll) close to the source. Often the nearest 
offset is 100 meters. For close measurements, more robust 
sensors may be needed and detector 'character' matching should 

20 be an important preliminary step. 

In US Patent No. 4,890,264 a method is proposed for suppressing 
non-uniformly distributed noise generated by surface wave 
propagation, wind, and machinery. A number of horizontally 

25 sensitive geophones are distributed amongst the conventional 
vertically oriented geophones. The outputs of the surface wave 
detectors are utilized in conjunction with an adaptive filter to 
cancel the effects of the surface wave interference. This method 
for the suppression of ground roll is inherently a multi- 

30 component method, and cannot be used in conjunction with single 
component acquisition. In addition, it neglects the fact that 
some seismic body wave energy also is detected by horizontally 
sensitive geophones, and this may cause signal cancellation. 
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In UK Patent Application GB-A-2273358 linearly constrained 
adaptive beamforming and adaptive interference canceling 
beamforming is used for ground roll suppression. In linearly 
5 constrained adaptive beamforming, signals measured by an array 
of geophones are filtered and summed so as to preserve signals 
incident from a preferred direction while suppressing 
interferences incident from other directions. In applying 
adaptive interference canceling, the moveout differential 

10 between the seismic reflections and the ground roll is used to 
form primary and reference channels. The filtering is performed 
using a continuously adaptive method such as the LMS ( least - 
mean-square) algorithm. The suggested application is in seismic 
while drilling, where the horizontal offset range is very small, 

15 so that the seismic reflections have an almost vertical angle of 
incidence and there is effectively a lot of data available from 
each source-receiver position since the roller cone drill bit 
used as the seismic source moves very slowly. The statistics of 
the noise then change very slowly, allowing stochastic gradient 

20 type of algorithms like the LMS to converge. However, in surface 
seismic experiments the ground roll present is often non- 
stationary and inhomogeneous . Therefore, stochastic gradient 
type of algorithms such as LMS may be too slow in converging. 

25 US Patent No. 5,237,538 proposes a method for removing coherent 
noise from seismic data. In this method, first the moveout 
characteristics of the noise are identified. Next, a space-time 
gate containing the noise defined and extracted, and the moveout 
removed to flatten the noise train. Amplitudes and time 

30 variations are then removed from the gate. The coherent noise is 
estimated using conventional stacking. A single -channel Wiener 
filter is. used to match the noise estimate to the noise in 
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the data trace containing signal -plus -noise . Having subtracted 
the filtered noise estimate, inverse amplitude scalars are 
applied to undo the effect of amplitude equalization. The 
signal is then moveout restored into the original seismic 
record. This particular method for removing coherent noise from 
seismic data is an application of the well-known technique 
called Postbeamf ormer Interference Cancelling. It has some 
particular shortcomings for application for ground roll 
attenuation. First, the signal always leaks into the ground roll 
estimate, especially for shorter arrays. There is always a 
component of the signal present at the reference channel which 
is co-located in time with the signal in the primary channel. On 
the other hand, when the arrays are allowed to be longer, the 
dispersion present in the ground roll make it difficult to 
achieve effective beamsteering . 

In marine seismic surveys, an acoustic source generates waves 
which travel through the water and into the earth. These are 
then reflected or refracted by the sub- surface geological 
20 formations, travel back through the water and are recorded by 
long hydrophone arrays which are towed near the surface of the 
water behind a seismic vessel . The hydrophones are mounted in 
streamer cables, or streamers. There are usually 1-12 streamers 
towed which are each several kilometers long. The streamers are 
25 made up of sections which may typically be 100-200 meters long; 
each section consists of hydrophones inside an outer skin 
which may be filled with oil, foam, or a more solid substance. 
Stress -wires and spacers form the internal skeleton of the 
streamer. 

30 

While the streamers are being towed behind the vessel, self- 
noise is generated due to a variety of sources. The lurching of 
the vessel, especially in rough seas, causes vibrations in the 
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stress-wires which interact with the connectors and the oil- 
filled skin, generating bulge waves (or breathing waves) which 
propagate down the streamers. The pressure variations are 
detected by the hydrophones, adding and corrupting the detected 
5 seismic signals. As the streamer moves through the water, 

boundary layer turbulence causes pressure fluctuations at the 
outer skin wall, which are again coupled to the hydrophones. 

Bulge waves may also be caused by eddy shedding under elliptical 
water motion about the streamer caused by wave action. 
A method of applying adaptive signal processing to the 
attenuation of bulge waves is described US Patent No. 4,821,241. 
There it is proposed to co- locate stress sensors with the 
hydrophones in the streamer. The stress sensors are responsive 
to mechanical stresses applied to the cable, but are 
substantially unresponsive to acoustic waves propagating in 
fluid media. The signal outputs from the stress sensors are 
combined with the signal outputs from the corresponding co- 
located hydrophones to cancel spurious signals due to bulge 
waves . 

Another method of applying adaptive signal processing to the 
attenuation of bulge waves was described is described US Patent 
No. 5,251,183. In this patent it is proposed to use an 
25 accelerometer secured between the lead-in section of the 
streamer and the hydrophone. Intra-shot and inter-shot 
accelerometer and hydrophone signals are recorded. The method 
utilizes inter- shot and intra-shot adaptive processing loops. 
The inter-shot adaptive processing loop derives inter-shot 
30 complex weights from inter- shot accelerometer signals and inter- 
shot hydrophone signals. The intra-shot adaptive processing loop 
models bulge wave noise in the intra-shot hydrophone signals by 
combining the inter-shot complex weights with intra-shot 
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accelerometer signals. Bulge wave noise attenuation is achieved 
by subtracting the intra- shot bulge wave noise model from the 
intra-shot seismic detector signals. 

5 SUMMARY OF THE INVENTION 

In accordance with the present invention, there is provided a 
method for filtering noise from discrete noisy seismic signals, 
comprising the steps of receiving signals using a plurality of 

10 receivers; determining a propagation characteristics of the 
signals with respect to receiver locations; and filtering 
received signals using an at least partially adaptive filter 
such that signals having a propagation characteristics other 
than the determined propagation characteristics are attenuated. 

15 The filtering step comprising the step of defining at least two 
independent sets of conditions (constraints) with a first set 
defining a desired (quiescent) response and a second set 
defining the propagation characteristics of signals to be 
preserved and the step of adapting filter coefficients of the 

20 filter subject to the independent sets of conditions 

(constraints) so as to minimize (optimize) the filter output for 
signals with a propagation characteristics other than the 
determined propagation characteristics. 

25 For the application of the invention, it is advantageous to 
define for the optimization process of the filter weights or 
coefficients a signal -dependent part (correlation matrix) and a 
signal -independent part. The signal -independent part usually 
comprises the constraints and is there often referred to as 

30 constraint matrix. Using this concept of a constraint matrix, an 
important aspect of the invention can be described as having 
within the constraint matrix a subspace which is defined by the 
desired quiescent response and one subspace which defines the 
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regions of the protected signal. By making these two subspaces 
orthogonal, filter weights can be found which simultaneously 
impose both restrictions upon the filter response. As the 
constraint matrix effectively reduces the degrees of freedom of 
5 the filter available for the adaptation process, this aspect of 
the invention can be described as splitting the total number of 
degrees of freedom into a first part which is available for the 
adaptation process and a second part which is used to define the 
constraints. The degrees of freedom assigned to the constraints 
10 are split among those constraints which defines the desired 

response and a second set defining the temporal and/or spatial 
spectral content or the propagation characteristics of the 
signals to be preserved. 

15 It is an advantage of the method to be not confined to narrow- 
band signals, but also applicable to wide-band seismic signals 
resulting in a filter that changes its response with the 
frequency of the input signal (s) . 

20 It is an important aspect of the invention that having derived a 
method of separating the desired quiescent response and the 
constraints relating to the region into orthogonal subspaces, 
any known method can be used for the adaptation process. Such 
known adaptive methods are known and described in the 

25 literature, e.g. LMS, RLS, LSL, FTF, etc.. 

According to a preferred embodiment of the invention, a filter 
bank comprising temporally and spatially local filters is used 
as the adaptive filter. 

30 

A filter bank can be defined as comprising M local multichannel 
adaptive filters with K channels, each of a length L. For most 
applications, the number L of coefficients is equal to or larger 
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than three. The number of channels K and of individual filters M 
are preferably two or more. 



The use of a filter bank for noise attenuation of seismic 
5 signals has been described in the International Patent 

Application W097/25632. However, the present invention does not 
require defining a reference channel in order to calculate the 
adapted filter bank coefficients. No noise estimate enters the 
adaptation process. Therefore, the present method can be applied 
10 to noise contaminated seismic signals, where there is no 

independent measurement or estimate of the noise available. 

According to one aspect of the invention, the coefficients of 
the filter are constrained such that its response corresponds to 
15 that of a beamformer with a specified look-direction. 

According to another aspect of the invention the constraints are 
set such that the filter preserves signals from a range of look 
directions or of defined regions of the f requency-wavenumber 
20 domain. The region can be pre-selected depending on the nature, 
more specifically on the apparent velocity of the seismic 
signals. Certain limits of the velocity, such as 1500 m/s, 
define regions in the f requency-wavenumber domain. 

25 A further aspect of the invention comprises the minimization of 
a cost functional using the approximation that the sum, weighted 
by window functions, of the output of adjacent filters of the M 
filters is equal when applied to the same signal in time regions 
where said window functions overlap. Preferably the method 

30 includes the step of multiplying M filtered estimates with 
temporal window functions. The application of the temporal 
window functions, and hence the resulting temporal windows, to 
the combined components ensures that the filtering process is 
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local in time and allows the method adapt ively to remove noise 
from the seismic data in accordance with a global optimization 
criterion. The data selection temporal window functions are 
preferably determined by two requirements, wherein the first 
5 requirement is that the sum over all windows at any given time 
equals unity and the second requirement is that only adjoining 
windows overlap. These requirements ensure that the global 
optimization of the filtered signal can be solved by use of an 
approximation in which for the sum over all time and all filters 
10 and all neighboring filters, the error function of a 
neighboring filter is replaced with the error function 
associated with the filter itself. 

The application of the data selection temporal windows decouples 
15 the equation required to solve the optimization of the filtered 
signal . 

According to yet another aspect of the invention, the response 
of the filter can be controlled by using a regularization 

20 parameter. The parameter as applied herein determines the 

relative weight of two components of the cost functional. One of 
the component of the cost functional can be defined as output 
power, while the other can be characterized as being essentially 
the white noise gain of the filter bank, i.e., the output of the 

25 filter in response to an input uncorrelated in time and space. 

The noisy signal may be pre-processed before being passed to the 
adaptive filtering means by dividing the signal into frequency 
bands using a reconstructing filter, for example a quadrature 
30 mirror filter. This allows a reduction in the number of data 
points to be processed and also allows a reduction in the number 
of coefficients in the adaptive filtering means as effectively 
reducing the bandwidth of the original signal. 
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The invention is applicable for two-dimensional (2D) and three- 
dimensional (3D) seismic surveys, and can be used in land 
seismic, marine seismic including sea bottom seismic, and 
transitional zone seismic. 

The method can be performed on stored data or on raw seismic 
data as it is acquired. Thus raw seismic data may be filtered 
according to the method at the data acquisition site. This 
ensures that a "cleaned" signal is available from the data 
acquisition site and may be downloaded directly from the site in 
this form. This reduces the amount of data sent for analysis 
off-site and reduces the costs and storage problems associated 
with accumulating sufficient quantities of noisy data for 
analysis off -site. The method can be advantageously applied to 
single-sensor recordings, i.e. to recordings prior to any group 
forming which combines the signals of two or more seismic 
sensors. 

20 Although the description of the present invention is based on 

seismic signal processing, it can be applied to sonic signals as 
used for example for well logging applications. Specific seismic 
applications include swell noise or streamer noise attenuation, 
including streamer noise attenuation in a cross -flow 

25 acquisition, attenuation of ground roll or mud roll or other 
coherent noise from marine, land, or transition zone data, 
seismic interference canceling, i.e. filtering noise using the 
full aperture of a mult i -streamer array, which is either towed 
in the water or deployed at the sea-bottom , or removal of sea- 

30 floor reflections from the notional source signature estimation, 
a technique described for example in the European Patent 
Application EP-A-066423. Other applications include noise 
suppression for various borehole seismic exploration methods, 
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where either noise or signal has preferential directions. Known 
borehole seismic methods include seismic-while-drilling (SWD) , 
vertical seismic profiling (VSP) , look-ahead and look-around 
while drilling. 



These and other features of the invention, preferred embodiments 
and variants thereof, possible applications and advantages will 
become appreciated and understood by those skilled in the art 
from the following detailed description and drawings. 



DRAWINGS 



FIG. 1 illustrates general elements of a seismic land 
15 acquisition; 



FIG. 2 shows a general block diagram of an adaptive 

beamformer in accordance an example of the present 
invention; 



FIG. 3 shows the example of a region preserved by an adaptive 

beamformer in accordance with an example of the 
pre sent invent ion ; 



25 FIG. 4 shows steps of defining a constraint matrix in 

accordance with an example of the present invention; 



FIG. 5 illustrates steps of using the constaint matrix in a 

process of filtering seismic recordings in accordance 
30 with the present invention; 



FIG 



6 . shows the example of a preserved region including a 
noise region suppressed by an adaptive beamforming 
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process in accordance with an example of the present 
invention; 



FIGs. 7A-C show examples of defining a preserved region or 

5 regions in accordance with the present invention. 

MODE(S) FOR CARRYING OUT THE INVENTION 

10 In the following the basic concepts underlying the invention are 
developed in detail. 

A typical land seismic acquisition is illustrated in FIG 1. 

15 A source 10 is activated, thus generating seismic waves 11, i.e 
acoustic waves with frequencies of less than 500 Hz. The waves 
travel though the interior of the Earth 12 and are reflected at 
various locations. Even though only one reflector 13 is shown, 
typically there are many reflectors, each reflecting a fraction 

20 of the seismic wave back to the surface. At the surface, seismic 
waves are recorded by seismic sensors 14 (geophones) . These 
sensors are spread along a line or in a two-dimensional pattern. 

As an example for a major source of noise, the travel path 15 of 
25 the so-called "ground roll" is shown. The ground roll is direct 
wave energy which propagates within layers close to the surface. 
It is one distinguishing feature of the ground roll to have 
different propagation characteristics than signals reflected 
from a deeper reflector layer: The ground roll reaches the 
30 sensors of the depicted sensor line one after the other. In 

contrast, the reflected seismic signal 16 from a reflector at 
great depth reaches all the sensors 14 almost simultaneously. 
When translated into the frequency- wavenumber domain, known as 
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f-k domain, the desired seismic signal therefore usually lies 
within a narrow cone around the f-axis (equivalent to small 
values of k) whereas the ground roll tends to have larger values 
of k. 

5 

Referring now to FIG. 2, there is shown a general block diagram 
of an adaptive beamformer in accordance with the present 
invention. It is assumed the presence of K sensors located at r k 
with k = 1,...,K. Each sensor k records a signal g k (n) with n = 

10 1,...,N. The letter 'n' is used as an index on discrete time 
samples. The sampling interval is Dt. The signals g k (n) are 
beamsteered using delays t k towards a general "signal direction" . 
This is the general direction from which the seismic signals are 
expected to arrive. The beamsteered channels x k (n) are processed 

15 by local multichannel adaptive filters to produce the output 
signal : 

M K L 2 

[i] y(n) = S Z Z h i (n)w ikv x k (n - v) 

i = lk = lv = -Li 

20 where w ikn (t) are the adjustable coefficients of the adaptive 
filters, hi(n) are the windows applied at the output end, M is 
the number of local multichannel adaptive filters (or the number 
of output windows) , and L = L x +L 2 +1 is the number of coefficients 
per channel. Here and below, a bar under a letter denotes a 

25 vector (small letter) or a matrix (capital letter) . 

Equation [1] can be rewritten as a (windowed) sum over a scalar 
product using a tap- input vector x(n) at time t defined as : 
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x(n) ■ [x x (n + L x ) , . . . , x x (n - L 2 ) , 
[2] x 2 (n + L x ) , . . . , x 2 (n - L 2 ) , 

x K (n + L x ) , . . . , x K (n - L 2 ) ] T 



and a tap-weight vector w A defined as 



Wj. = [w il( _ Li)/ . . • , W ilL2 , W i2( _ Li)7 . . . , W i2L2 , 



5 [3] _ T 



W 



iK(-Li)' " ' ' W iKL 2 J 



Using definitions [2] and [3] , equation [1] becomes 



M M 

[4] y(n) = 2 hitnJWiXjJn) = £ h^rOx (n)Wi . 

i = l . i = l 

10 

Equations [1] and [4] describe how to find the beamformer or 
filter bank output once the M tap-weight vectors w A have been 
specified. These vectors are computed as the solution of an 
optimization problem which is described below. 

15 

The optimization problem is defined as 



8 2 

[5] min J = min \J 1 + — J 2 

w.i,...,w M 3£i,...,w M 



KL 



20 subject to constraints 



[6] C T Wi = f 



where i = 1,2, . . . , M and 
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[7] 




Z Y (n) and 



n = l 



[8] 




Z||wJ 2 £ hi (n)||x(n)|| 2 , 



i=l n=l 



5 

KL is the total number of filter coefficients, and | | . | | denotes 
the L 2 norm. This cost functional is a linear combination of the 
output power of the beamformer (the first term in eq. [5]), and 
the so-called "white-noise gain" of the beamformer weighted by 

10 the input power (the second term in eq. [5]), The relative 

weight of the two terms is adjusted by the d 2 term. Including the 
"white-noise gain" of the beamformer in the cost functional is 
intended to increase the beamformer robustness in the presence 
of signal modeling uncertainties (sometimes called 

15 perturbations,) and numerical correlation between the signal and 
the noise. 

Equation [6] describes Q linear constraints on the admissible 
solutions to the optimization problem. Here, the KLxQ matrix C 
20 is the constraint matrix, and the Q-vector f. is the response 
vector. The actual design of the linear constraints are 
discussed below. 

A possible solution of the optimization depends on imposing the 
25 following two constraints on the window functions h t (n) : 



[9] 



Z ^(n) = 1 



i = l 



for n 



2 



N, and 
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[io] h^rOhjdi) = 0 

for j <> i-l,i,i+l. The first constraint ensures that the filter 
5 bank is equivalent to the single filter case if all the local 
filters (Wi)are identical. The second constraint ensures that 
the windows have compact support. 

The optimization problem can be to a large extend decoupled 
10 using the second condition(eq. [10]), and the approximation 

I X S hi(n) hj(n) w? x(n) x T (n) Wj « 

n i j=i-l,i+l 

[11] 

SI Z hi(n) hj(n) x(n) x T (n) w ± 

n i j = i-l,i + l 

The approximation of equation [11] requires that neighboring 
15 filters produce similar results when applied to the same input 
data in time regions where adjacent windows overlap, instead of 
requiring that neighboring filters are similar on a point -by- 
point basis. Thus, the approximation is similar to requiring 
that the integral of two functions are close, rather than the 
20 functions themselves. 

With this approximation, the first term of the cost functional, 
Jj , becomes 

M 

25 [12] Ji = Z wjOiWi 

i = l 



with 
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[13] 



Oi = X hi(n)x(n)x T (n) . 



n 



The second term in the cost functional can be rewritten as 



[14] 



M 
2 ~ 

i = l 



2 tr\z h^njxfajx^n) 

n = l 



with "tr" denoting the trace of a matrix. 



10 Combining Equations (5) , (12) , (14) , and reorganizing the terms, 
the total cost functional can be written as: 



[15] 



M 5 

J = Z + — tr(Oi)I 



Zu Wil 

i = l I 



KL 



15 where I denotes the KLxKL identity matrix. The decoupled 
optimization problem can be solved for each of the M time 
windows subject to the constraints [6] . Using the method of 
Lagrange multipliers, the optimal tap-weight in each window is 



given by 



20 



[16] 



Wi* = ^.i' 1 C(C T O^ 1 0' 1 f , 



with 



25 [17] 



5>i = + 



KL 



tr(Oi)I, 
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The second term of the modified local correlation matrix F t " can 
be regarded as a regularization term with d 2 as the 
regularization parameter. In array signal processing literature, 
regularization of correlation matrices with the addition of a 
5 scaled identity matrix has been suggested to increase robustness 
in the presence of perturbations, in the context of narrow-band 
beamforming. Here, the cost functional [5] includes the 
regularization term from the beginning leading to a 
generalization for wide -band adaptive beamforming. Hence, the 
10 filter response changes as a function of the frequency of the 
signal . 

When the input data to the beamformer is characterized by 
spatially and temporally uncorrelated (or white) noise, both the 
15 correlation matrix F t and the modified correlation matrix F t ~ 
become proportional to the identity matrix. In this case, the 
optimal weight vector becomes 



[18] Wi* a w q = C(C T C) 1 f . 

20 

The weight vector Wg is called the quiescent solution to the 
optimal beamformer problem, and the corresponding response is 
known as the quiescent response. Note that the quiescent 
solution depends entirely on the constraint matrix C and the 
25 response vector f. . 



The optimal weight vector w A approaches the quiescent weight 
vector even for general noise fields as the regularization 
parameter d 2 is increased. In this case, the modified correlation 
30 matrix Z± approaches the identity matrix (cf. [17]) . The 

■ 

regularization parameter d 2 therefore weights the optimal 
solution between a solution that is entirely dependent on the 
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received data, and a solution that is independent of the data. 
For d 2 = 1, both solutions are equally weighted in the sense that 
their corresponding correlation matrices have equal trace value. 
In situations where the perturbations are higher, i.e. the 
5 assumptions about the seismic acquisition geometry do not 

exactly hold, finding a beamformer response with a higher level 
of regularization can give more robust results. 

Another aspect of the invention relates to the design of linear 
10 constraints (eq. [6]) to be imposed on the beamformer. 

One type of linear constraints that can be imposed on the 
beamformer are those designed to preserve seismic signals 
incident from a target direction, while suppressing 

15 interferences incident from other directions. Steering delays t k 
as those shown in FIG. 2 define a single "look-direction" . 
Signals incident in this direction are in phase, and for these 
signals the system can be considered as a single FIR (finite 
impulse response) filter. The values of the coefficients for 

20 this equivalent processor are equal to the sums of the 

corresponding coefficients in the adaptive processor. Each local 
beamformer w t consists of the adaptive filters w n , w i2 , . . . ,w iK 
processing data from each channel, and a summing unit. The sum 
of the individual filters w ix , # i2 ,...,w iK is constrained to give 

25 w^, which is the desired response for signals incident in the 
look-direction, e.g., a unit pulse in the look direction. The 
quiescent response then becomes that of a fixed- weight 
beamformer with single equal weights for all elements. In the 
frequency- wavenumber domain, this corresponds to a sync function 

30 that is constant in the f direction. Therefore, for increasing 
values of the regularization parameter d 2 , the beamformer 
preserves, signals incident not only from the look direction, but 
also from neighboring directions. 
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As discussed in the last section, using single look-direction 
constraints and regularization, it is possible to preserve 
signals incident from directions near the look direction. While 
5 this approach is useful and sufficient for many applications, it 
is desirable to derive more general linear constraints that will 
satisfy the requirements in any seismic data acquisition 
situation more directly. 

In narrow-band beamforming, different generalized constraint 
design approaches are known. Derivative constraints are used to 
influence response over a region of the response space 
by forcing the derivatives of the beamformer response at certain 
points of the response space to be zero. Eigenvector constraints 
are based on a least squares approximation to the desired 
response, and are usually used to control the beamformer 
response over regions of the response space. Generalization of 
these methods to wide -band beamforming problems have shown that 
while they provide a good response in selected regions of the 
response space, they can generate unacceptably high sidelobes in 
other regions. 

For the present invention, the requirements of the generalized 
constraint design are to impose an arbitrary quiescent response 
25 on the beamformer and to make sure that certain areas in the 
frequency-wavenumber domain are entirely controlled by the 
quiescent response. These requirements have been established 
with the following functional objectives in mind: 

30 - accommodation of an arbitrary range of apparent signal 
velocities; 

- increased robustness to perturbations; 

- capability to use larger arrays; 



10 



15 



20 
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- being able to run the adaptive beamformer with a lower 
regularization level (d 2 ) , hence achieving higher noise 
attenuation; and 

- achieving higher noise attenuation for a given level of 
5 regularization by the appropriate design of the quiescent 

response . 

To impose an arbitrary quiescent response on the beamf ormer, use 
can be made of the fact that the linear constraints [6] define a 
10 Q-dimensional hyperplane in a KL- dimensional space. Equation 
[18] shows that the quiescent weight vector Wg is the minimum 
norm solution to Equation [6], i.e., it is the shortest vector 
from the origin to the hyperplane. 

15 Equation [18] also shows that Wg is a member of the subspace 
spanned by the columns of the constraint matrix C. The columns 
of C are in general independent (otherwise some constraints 
would be redundant) , thus they can be chosen to be orthogonal 
without loss of generality. After defining a desired quiescent 

20 weight vector w^, this suggests the following forms 
for the constraint matrix C and the response vector f : 

[19] C = [aw qd , d] 

25 and 

[20] f = 

with the condition 

[21] p = a 




30 
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where D is a KLx(Q-l) matrix whose columns are orthogonal to w^. 
The exact form of the matrix £ is described below. With C and f 
chosen according to [19] and [20] , respectively, it can be shown 
that the desired weight vector equals the quiescent response 
5 vector (eq. [18]). 

After defining the first column of the constraint matrix C and 
the response vector £ to impose the quiescent weight vector, the 
definition of the matrix 2 which is a part of C is derived in 
10 the following sections. 

In a seismic acquisition, reflection signals that should be 
preserved can be considered as a linear combination of plane 
waves with associated frequency and wavenumber values from a 

15 known region of the f requency-wavenumber space. This region, 
which is denoted A in FIG. 3, depends on the particular 
acquisition geometry, but is usually a cone around the frequency 
axis. One possible example of a preserved region in the 
frequency -wavenumber domain is shown in FIG . 3, where A is 

20 chosen so as to include all signals of apparent velocity of +/- 
1500 m/s or more. In the present example, the beamformer 
response in region A should be controlled entirely by a 
quiescent response which preserves the signal. 

25 The set S A of seismic signals to be preserved by the filtering 
process is defined by 

[22] S A = {s(t, r) : s(t, r) = JJ A df dk S(f , k^ 2 *^"^} 

30 as composites of plane waves with associated frequencies and 
wavenumber values from a region A, where S(f,k) is the complex 
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Fourier amplitude corresponding to the plane wave component of a 
signal with frequency "f and wavenumber ]£. 



Using [22] , the tap- input vector [2] can be rewritten as 



[23] 



s(n) = JJ (f #k€AJ df dk S(f , k)e j2nfnAt d(f , k) , 



10 



15 



with £(f ,k) being defined as the array steering vector 
corresponding to the plane wave component specified by 
particular frequency f and wavenumber )£. It is noteworthy that 
in contrast to the example described above no time delays t have 
been introduced into the signal path to steer the filter 
response. The array steering vector can be written as a 



Kronecker product : 



-j27t]£'£i 
-j27ik-£ 2 



j27ifL 1 At 



[24] 



d(f, k) = 



-j2rck-r K 



-j27tfL2At 



Using [4] , the response of the beamformer to the signal tap- 
input vector s(n) is 

20 

[25] y(n) = 2 hi(n)Jj A df dk S(f , k) T e j27lfnAt d T (f , k) Hi ■ 

i = l 

For the beamformer response to be the same for both the optimal 
weight vector ^ and the quiescent weight vector w^, and further 
25 requiring that this equality to hold for all signals s(t; r) of 
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the preserved region, i.e. signals with arbitrary associated 
Fourier coefficients S(f, k) such that (f, k) is in A. This 
requires that 

5 [26] d T (f , £>Wi* = d T (f , w q , V(f , k) e A . 

By decomposing the optimal weight vector into a fixed weight 
portion equal to the quiescent weight vector and an adaptive 
weight portion according to a solution known as "generalized 
10 sidelobe canceller" (GSC) , it can be shown that the last equation 
is equivalent to requiring that d(f ,k) lies in the column space 
of the constraint matrix £. 

It is therefore a further object of the present invention to 
15 find a efficient, i.e. preferably low rank, basis for the space 
of steering vectors d(f ,k) . However, a scalar multiple of w^ has 
already been installed as the first column of C, we actually 
need to find a low rank basis for the part of this space that 
lies in the orthogonal complement subspace of Wq d . The projection 
20 of £(f ,k) onto the orthogonal complement of w^ is the projected 
steering vector: 

[27] 3(f, k) s d - P Hqd ) d(f, k) , 

25 where the expression in parentheses is the orthogonal complement 
projector with respect to Wq d with 

■ 

t281 £w ad = ^qdteqd^qd)" 1 ^/ 
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Using the. fact that any KL- dimensional d"(f,k) can be written as 
a linear combination of orthonormal vectors {v x , . ., v^} , 
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KL 



[29] 



3(f/k) = Icr D (f,k)v p - vg(f,k) , 



P =i 



a rank P (P < KL ) approximation of the projected steering 
5 vectors is obtained by 



[30] 



3 P (f , k) = I 3 D (f , k)vp a Va(f , k) , 



P =i 



where 



10 [31] 



a(f , Is) = [ajf , Is) / • • • / cr p (f , k) , 0, . . . , o] T . 



To derive an efficient rank P representation of d/(f,k) for any 
(f ,k) in region A, an error functional with respect to the L 2 
norm is defined as 



15 



[32] 



Hp - JJ. df dkS(f,k) - S P (f , k) 



Using the correlation of all projected steering vectors in 
region A of the f requency-wavenumber space given by 



20 



[33] 



R A = JJ. df dk 3(f , k)2 H (f , k) 



the error functional can be expressed as 



KL 



25 [34] 



p=p + l 
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The superscript W H" denotes the conjugate transpose of a vector 
or matrix. 

The optimum set of ordered basis vectors v can be found by 
minimizing the cost functional m p subject to the constraint that 
Vp H v p = 1, with 1 <= p <= KL. Using Lagrange multipliers, the task 
is to minimize 

KL r 1 

[35] Z Y«R A v p - Vv£v p - DJ ■ 

p=p+l 

By taking the gradient with respect to Vp and setting it to zero, 
the optimal basis vectors v^...,^ are found as the eigenvectors 
of R^fwith respect to the eigenvalues l p ) - The missing part D of 
the constraints matrix £ (cf. [19]) can now be defined as the 
principal eigenvectors of R^: 

[36] D = [V-l, . . . , V p ] . 

Note that the steering vectors d(f ; 3s) are in general complex 
valued. Therefore, their correlation matrix E*" over a general 
region A in the f requency-wavenumber space is complex valued, 
making the eigenvectors of R^' , hence the columns of C also 
complex valued. However, in seismics the signals are real valued 
signals which have complex conjugate Fourier coefficients. 
Therefore the types of A regions that are of interest are always 
symmetric in the f requency-wavenumber space with respect to the 
origin. The resulting matrices (R^, C)are then all real valued. 

The above described expansion of the projected steering vectors 
£T(f; 3i) is analogous to the Karhunen-Ldeve expansion. While the 
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original Karhunen-Ldeve expansion is for a random vector, the 
expansion presented here is for a deterministic set 
of vectors- This is reflected in the way the approximation error 
functional m p is defined, cf . [32] . 

5 

The covariance matrix of steering vectors, similar to the 
correlation matrices defined in [33] was first introduced in by 
K.M. Buckley, IEEE Trans. Acoust . Speech Signal Processing, Vol 
ASSP-35, 249-266, March 1987, but was then heuristically 
10 defined within a stochastic framework, assuming zero mean 
signals and using a narrow-band representation of wideband 
signals. In the description of the present invention, the 
correlation matrix has been derived from first principles within 
a deterministic framework. 

15 

The main steps of the generalized constraint design method are 
shown in the flow diagram of FIG 4 . They include : 

- specification of the desired quiescent weight vector, w^, 
which defines the first column of the constraint matrix; 

20 - specification of the signal protection region A in the 
frequency- wavenumber space; 

- computation of R^' , the correlation matrix of all the projected 
steering vectors in region A; and 

- determination of the principal eigenvectors {v l# , . ., v^} of 
25 as the remaining columns of the constraint matrix. 

Having computed these, the constraint matrix is specified as 




and the response vector as 
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[38] f = 

using a convenient choice for a and b in eq. [21] 

5 The specification of the desired quiescent weight vector to form 
a desired quiescent response is essentially a non-adaptive 
multidimensional filter design problem, for which 
many techniques exist. Reference can be made for example to 
handbooks such as W. Chen (ed.), "The Circuits and Filters 
10 Handbook", IEEE and CRC Press, 2732-2761 (1995), D.E. Dudgeon 

and R.H. Martinez, "Multidimensional Digital Signal Processing", 
Prentice Hall (1984), or J.S. Lim, Two-Dimensional Signal and 
Image Processing, Prentice Hall (1990) . 

15 Once the constraint matrix is defined, digitized recordings of 
single seismic sensors can be filtered to generate a single 
sensor recording with reduced noise. The "cleaned" recording can 
then be used in further processing steps, such as group forming, 
stacking, velocity analysis, moveout correction etc., known in 

20 the art to ultimately generate a representation of subterranean 
formations. These steps are outlined in the flow diagram of FIG. 
5. As however the details of these steps (with the exception of 
the filtering step) are of no particular concern regarding the 
present invention, are detailed description thereof is omitted 

25 herein. 

The following section presents alternative methods of 
efficiently defining the region A (see FIG. 3) protected by the 
quiescent response. 

30 
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For seismic applications, the quiescent weight vector could be 
designed such that over the region A in the f requency-wavenumber 
space the response is near unity, thus preserving the seismic 
signals in that region. In regions of the f requency-wavenumber 
5 space where the noise is expected to be present, the quiescent 
response should have low values, so that even when 
regularization is used, high performance can be achieved. 

The constraint design process can be extended as described in 
10 the next section. The constraint design outlined above resulted 
in a low rank basis of projected steering vectors in A. The 
objective was to preserve all signals in the preserved region A 
without any reference to their relative strength. This is 
reflected in the choice of the error functional m p defined in 
15 [32] . In many applications this choice makes sense, where it is 
desirable to protect signal components which have much lower 
amplitude then other signal components. On the other hand, in 
some other applications it may be desirable to minimize the 
power of the overall signal distortion. 



This extension of the above described method is reflected by a 
generalization of the definition of m p to 



20 



[39] 



Hp s If A df ^ S(f , k) |3(f , k) - 3 P (f , k)|| 2 . 



25 



The correlation matrix of the original steering vectors in A 



then becomes 



[40] 




30 
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A correlation matrix of the projected steering vectors can be 
derived from [40] using the orthogonal complement projector 
(eqs. [27] , [28] ) : 

5 [4i] £ A = (I - Pw qd )R A G - Ew qd ) ■ 

In some applications, it may be desirable to reserve additional 
regions of the f requency-wavenumber space over which the 
beamformer response would be controlled entirely by the 

10 quiescent response. In FIG. 6, there is shown an example of 

region A in the f requency-wavenumber space which includes signal 
protection region Al and noise region A2 . If it is known that in 
a certain environment there is almost always a coherent noise 
component occupying the region A2 , then it may be beneficial to 

15 design the quiescent weight vector Wq d which would put deep nulls 
in that region, i.e, which suppresses any signal from a 
direction (f,k) where (f,k) is element of A2 . Then, the adaptive 
weights would concentrate on attenuating the noise in the 
remaining regions of the f requency-wavenumber space. 

20 

Employing the generalized constraint design methodology, the 
quiescent response of the beamformer can be specified exactly, 
and the beamformer response over a pre -specified region A in the 
f-k space can be constrained to approximate the quiescent 

25 response to an arbitrary extent. The accuracy of this 

approximation is controlled by the user parameter P, which is 
the number of principal components of R*~ • As P is increased, 
more and more degrees of freedom of the beamformer are fixed, 
and the adaptive degrees of freedom are reduced. As P approaches 

30 KL-1 the beamformer response approaches the quiescent response 
not only in A, but over the entire f-k space regardless of the 
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noise field characteristics or the regularization parameter 
used. 

When the region A is relatively large, the distribution of the 
5 eigenvalues of R^' may be only slowly decreasing, so that the 
number of principal components required to adequately represent 
the space of steering vectors in A may be high. If only a small 
number of principal components are used in order to keep more 
degrees of freedom adaptive, the response of the beamformer may 
10 deviate from the quiescent response in A significantly. This may 
happen even over the original look direction indicated by k = 0. 
This is in contrast with the original single "look-direction" 
constraints described before, which guarantee that the response 
over the k = 0 line is identically unity for all frequencies. 

15 

In order to guarantee that the response of the beamformer is 
exactly that of the quiescent response over k = 0 as in the case 
of single look-direction constraints while constraining the 
response arbitrarily tightly in the rest of region A, the region 
20 A is preferably partitioned and each section is treated 
separately. 

Methods of pursuing this approach are shown in FIGs. 7A-C, where 
A is the region as shown in FIG 3. In the following examples, 
25 this region A is partioned into sections Al and A2 as 
illustrated in FIGs. 7A and 7B. 

It can be shown that if Al is the k = 0 line in the f-k space, 
then the space of projected steering vectors over Al can be 
30 spanned by L-l eigenvectors of R^" (see eq. [42] below) , as long 
as the quiescent response Wg has unity response over Al , as 
usually is the case. That is, the matrix 



- 34 



[42] g A1 = j A1 df ace, q) a H (f , o) . 

has rank L-l. This is an interesting observation, as it neatly 
ties the single look-direction constraints and the generalized 
5 constraints. If the signal protection region A is chosen as 
Al , the quiescent response is chosen as 



[43] W g = 



10 and the columns of matrix 2 in the constraint matrix C are 
chosen as the L-l eigenvectors of R^* corresponding to the 
non-zero eigenvalues, then the linear constraints set by the two 
methods become identical. 

15 Having above defined the region Al as the k axis, there are 
various possibilities for defining A2 . One would be to define 
A2 as the rest of the region A as in Figure 7A. Another way to 
define A2 is shown in Figure 7B. In this case A2 is basically 
the boundary of the original region A. 

20 

In general, if the number of array elements K is small and the 
original region A has a relatively short wavenumber extent, 
setting the additional constraints as in Figure 7B would be 
preferred. When the number of sensor elements is small., the 
25 variability of the beamformer response as a function of 
wavenumber is limited, and constraining the response 




10 
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more strictly at a few points for each frequency may be more 
effective than constraining it over all the points of a wide 
region. If the number of array elements is large or the region 
A is wide in wavenumbers, defining A2 as in Figure 7A would be 
preferred. It is of course possible to generalize even further 
to combine the two types of additional constraints as shown in 
Figure 7C giving rise to three sections denoted Al, A2 and A3 of 
region A. 

When the protection region A is partitioned using any of the 
approaches above, the constraints can be computed defining the 
response vector as in eq. [38] , where Q. has appropriate length, 
and the constraint matrix Q is given as 



15 [44] 



C = 



™qd / 



™qd 



/ SaI' 2a2 



Here, 2 Al is the matrix whose columns are the L-l principal 
eigenvectors of R^' (cf. eq. [42]) 



20 [45] 



Sai = G " E 



w gd )R A i(i-P 



) 



25 



using the projection matrix with respect to w^ defined in eqs 
[27] , [28] . ^ is the matrix whose columns are the principal 



eigenvectors of 



[46] 



Sa2 = (I- 2c a1 )Ra2<I " EgJ = 

a - p 8ai ) g - p Wgd ) r A2 (i - Pw gd ) a - E Dai 



) ' 



where 
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[47] C A1 = [^W qd / ||w qd || 2 ,D A1 j, 



[48] Pg A1 = £ A1 (C A1 C A1 ) C A1 



[49] Ed a1 = D A1 (P A1 T D A1 ) 1 D A1 T < and 



[50] R A2 = |L df dk d(f , Is) d H (f , k) . 



Equation [46] shows two alternative ways of computing R^". These 
10 equations have been written for partitioning of A into two 
subregions, but can be generalized to more sub-regions. 

By partitioning the original signal protection region A as 
described, choosing Al as the k=0 line, setting the quiescent 

15 response w as in [43] , and constructing the columns of matrix D 
in the constraint matrix C are chosen as the L-l eigenvectors of 
Bu" corresponding to the nonzero eigenvalues, any additional sub- 
regions A2 , A3 , etc. would give additional control over 
beamformer response with respect to the look- direct ion 

20 constraints. 

For some applications, it may be useful to reduce the degrees of 
freedom used by the adaptive beamformer. In the so-called 
partially adaptive beamformer, only a portion of the available 

25 degrees of freedom are used adapt ively. The main advantages of 
reducing the adaptive degrees of freedom are reduced 
computational cost and improved adaptive convergence rate. The 
primary disadvantage of partially adaptive beamforming is a 
degradation in the steady state interference cancellation 

30 capability of the beamformer. Therefore, the objective of 
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partially adaptive beamformer design is to reduce the number of 
adaptive weights without significantly degrading the performance 
of the adaptive array, 

5 Previous partially adaptive methods includes numerical 

techniques for approximately minimizing the average generalized 
sidelobe canceller (GSC) output power for a desired number of 
adaptive weights, where the average is over a range of jammer 
parameters The present invention uses a method which is based 

10 on a design method described by H. Yang and M. A. Ingram, IEEE 
Trans. On Antennas and Propagation, Vol.45, 843-850, May 1997. 
It also attempts to minimize the average GSC output power, but 
under a constraint that the reduced-dimensional solutions for 
all of the scenario trials lie in the same subspace. This 

15 constraint makes it possible to use a singular value 

decomposition to get the rank-reducing transformation, thereby 
simplifying the optimization problem. 

The generalized sidelobe canceller solution can be written as 
20 (cf. [18]): 

[51] = W q - B w ai , 

where B is a KLx(KL-Q) full rank matrix whose columns span the 
25 orthogonal subspace of the constraints matrix C and is known as 
the blocking matrix. The vector w^ is the KLxQ dimensional 
adaptive part of the optimal weight vector and is given by 

[52] w ai = (B^iBT 1 B T $iW q . 

30 

The partially adaptive GSC achieves a smaller number W of 
adaptive weights, through the use of a (KL-Q)x W linear 



- 38 - 

transformation T following B prior to adaptive weighting. The 
partially adaptive optimal weight vector can be expressed as 

[53] Wi* = W q - BTw pi , 

5 

where the W-dimensional adaptive part of the optimal weight is 

[54] w pi = (T T B T $ i BT)~ 1 T T B T $iW q . 

10 It is now the aim to choose T which minimizes the interference 
and noise output power over a set of likely interference 
scenarios- These scenarios can be characterized by different 
parameters like the number of interferers, interferer 
directions, interferer spectral densities, white noise levels, 

15 etc. The applied method can be summarized as follows: 

- for each random outcome qj from a distribution of scenario 
parameters, compute the full -rank optimal adaptive weight vector 
w.< from [52] and the transformed weight vector a given by 

20 

[55] a(ej) = USU T w ai (9j) , 

where 

25 [56] B T 3>i(9j)B = US 2 U T 

is the eigendecomposition of Eftqj); 



- store vectors Wai(qj)and a(qj) into the matrices W and A, 
30 respectively; 



- compute the singular value decomposition of & to get from 

[57] A = U h ¥i ■> 

and 

- derive £ as the first W columns of WA*U^ , where the superscript 

indicates the pseudoinverse . 

In most seismic surveys, the noise such as ground roll or swell 
noise occupies only a fraction of the temporal bandwidth 
available. For example in a land seismic survey, the Nyquist 
frequency is 2 50 Hz, while most of the ground roll energy is 
under 30 Hz. Concentrating filtering efforts to the frequency 
band where the noise resides is desirable to reduce 
computational cost . 

One means of achieving this aim involves adding QMF (quadrature 
mirror filter) perfect reconstruction filter banks, described 
for example by P.P. Vaidyanathan, in "Multirate Systems and 
Filter Banks, Prentice Hall, 1993 or by M.J.T. Smith and T.P. 
Barnwell III, in: IEEE Trans. Acoust. Speech Signal Processing, 
Vol. ASSP-34, 434-441 (1986) to the seismic noise and 
interference suppression system using adaptive multichannel 
filter banks, as described above. Two filter banks are used in 
this system. The QMF filter bank is used to decompose the traces 
into frequency bands and decimate before adaptive filtering is 
applied, and is subsequently used for resynthesizing the 
original signal. The multichannel adaptive filter bank is the 
heart of the system performing the actual filtering for noise 
suppression. Using the perfect reconstruction filter banks to 
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decimate reduces the number of points to be processed and also 
allows reduction in the number of coefficients in the adaptive 
filters, bringing in significant savings in CPU time and memory 
requirements. 
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CLAIMS 

1. Method for filtering noise from discrete noisy (seismic) 
signals, said method comprising the steps of 
5 - receiving signals using a plurality of receivers; 

- determining a propagation characteristics of said signals with 
respect to receiver locations; and 

- filtering received signals using an at least partially 
adaptive filter such that signals having a propagation 

10 characteristics other than said determined propagation 
characteristics are attenuated; 
said filtering step comprising the steps of 

- defining at least two independent sets of conditions with a 
first set defining a desired response and a second set defining 

15 said propagation characteristics of signals to be preserved ; 
and 

- adapting filter coefficients of said filter subject to said at 
least two independent sets of conditions so as to optimize the 
filter output for signals with a propagation characteristics 

20 other than said determined propagation characteristics. 

2 . The method of claim 1 wherein the propagation 
characteristics of a signal is determined by its temporal and 
spatial spectral content. 

25 

3. The method of claim 1 wherein the propagation 
characteristics of a signal is determined by its )g value and 
frequency. 

30 4. The method of claim 1 wherein the propagation 

characteristics of a signal is defined as a region or as regions 
in the f requency-wavenumber domain. 
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5. The method of claim 1 wherein the conditions define a 
constraint matrix comprising at least two mutually orthogonal 
subspaces with a first subspace defining the desired response of 
the filter and a second subspace defining the propagation 

5 characteristics of the signals to be preserved with said matrix 
being applied during the adaptation process of filter 
coefficients . 

6. The method of claim 1 wherein the conditions define a 
10 constraint matrix using a quiescent response vector and 

principle eigenvectors of a correlation matrix made of steering 
vectors defining at least one region of the frequency- wavenumber 
domain wherein the signals are preserved. 

15 7. The method of claim 1 defining the conditions such that 
signals with pre -determined propagation characteristics are 
suppressed by the filter. 

8. The method of claim 1, defining further sets of conditions 
20 with one set of conditions forcing the desired response of the 
filter to zero and one set of conditions defining the 
propagation characteristics of signals to be suppressed by the 

filter, 

25 9- The method of claim 1 wherein the filter comprises M 

temporally local filters, said filters forming a filter bank, 
and M being a number equal to or larger than two. 

10. The method of claim 9 including the step of multiplying M 
30 filtered estimates with temporal window functions (h^n)). 
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11. The method of claim 10 wherein the temporal window 
functions are characterized by the requirement that only 
adjoining windows overlap. 

5 12. The method of claim 1 wherein in the adaptation process a 
cost functional (J) representing the output power of the filter 
and its white noise gain is optimized. 

13. The method of claim 12 wherein the relative weight of the 
10 output power and the "white noise gain" in the cost functional 

(J) is adjustable (using a parameter d) . 

14. The method of claim 12 wherein the cost functional (J) is 
minimized using the approximation that the sum, weighted by 

15 window functions, of the output of adjacent filters of the M 

filters is equal when applied to the same signal in time regions 
where said window functions overlap. 

15. The method of claim 1 comprising the further step of 

20 dividing filter coefficients into a fixed part and a adaptive 
part . 

16. The method of claim 1 comprising the further step of 
splitting the signals into at least two frequency bands, 

25 filtering at least one of said at least two frequency band, and, 
before further processing, recombining said frequency bands to 
regain noise attenuated seismic signals. 

17. The method of claim 1 wherein the discrete noisy signals 
30 used as input are recordings from individual seismic sensors 

prior to any group forming techniques. 
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18. The method of claim 1 comprising the further steps of 
positioning at least one seismic source and a plurality of 
seismic receivers at appropriate positions in a land, marine or 
transitional zone environment; activating said at least one 

5 source to propagate energy through subterranean formations, and 
using said receivers to measure energy as single sensor 
recordings, converting said single sensor recordings into 
discrete seismic signals; and transmitting said signals as input 
to the filter. 

10 

19. A method of processing seismic signals in a to generate a 
representation of a subterranean formation wherein noise 
present in said signals is removed by filtering said signals 
using the steps of - receiving signals using a plurality of 

15 receivers; 

- determining a propagation characteristics of said signals with 
respect to receiver locations; and 

- filtering received signals using an at least partially 
adaptive filter such that signals having a propagation 

20 characteristics other than said determined propagation 
characteristics are attenuated; 
said filtering step comprising the steps of 

- defining at least two independent sets of conditions with a 
first set defining a desired response and a second set defining 

25 said propagation characteristics of signals to be preserved ; 
and 

- adapting filter coefficients of said filter subject to said at 
least two independent sets of conditions so as to optimize the 
filter output for signals with a propagation characteristics 

30 other than said determined propagation characteristics. 
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Amendments to the claims have been filed as follows 

1. Method for filtering noise from discrete noisy (seismic) 
signals, said method comprising the steps of 
5 - receiving signals using a plurality of receivers; 

- determining propagation characteristics of said signals with 
respect to receiver locations; and 

- filtering received signals using an at least partially 
adaptive filter such that signals having propagation 

10 characteristics other than said determined propagation 
characteristics are attenuated; 
said filtering step comprising the steps of 

- defining at least two independent sets of conditions with a 
first set defining a desired response and a second set defining 

15 said propagation characteristics of signals to be preserved ; 
and 

- adapting filter coefficients of said filter subject to said at 
least two independent sets of conditions so as to optimize the 
filter output for signals with propagation characteristics other 

20 than said determined propagation characteristics. 

2 . The method of claim 1 wherein the propagation 
characteristics of a signal is determined by its temporal and 
spatial spectral content, 

25 

3 . The method of claim 1 wherein the propagation 
characteristics of a signal is determined by its k value and 
frequency. 

30 4. The method of claim 1 wherein the propagation 

characteristics of a signal is defined as a region or as regions 
in the f requency-wavenumber domain. 

5. The method of claim 1 wherein the conditions define a 
35 constraint matrix comprising at least two mutually orthogonal 



5. The method of claim 1 wherein the conditions define a 
constraint matrix comprising at least two mutually orthogonal 
subspaces with a first subspace defining the desired response of 
the filter and a second subspace defining the propagation 

5 characteristics of the signals to be preserved with said matrix 
being applied during the adaptation process of filter 
coefficients . 

6, The method of claim 1 wherein the conditions define a 
10 constraint matrix using a quiescent response vector and 

principle eigenvectors of a correlation matrix made of steering 
vectors defining at least one region of the f requency-wavenumber 
domain wherein the signals are preserved. 

15 7. The method of claim 1 defining the conditions suck that 
signals with pre -determined propagation characteristics are 
suppressed by the filter. 

8. The method of claim 1, defining further sets of conditions 
20 with one set of conditions forcing the desired response of the 
filter to zero and one set of conditions defining the 
propagation characteristics of signals to be suppressed by the 
filter. 

25 9. The method of claim 1 wherein the filter comprises M 

temporally local filters, said filters forming a filter bank, 
and M being a number equal to or larger than two. 

10. The method of claim 9 including the step of multiplying M 
30 filtered estimates with temporal window functions (h^n)). 



11. The method of claim 10 wherein the temporal window 
functions are characterized by the requirement that only 
adjoining windows overlap. 

12. The method of claim 1 wherein in the adaptation process a 
cost functional (J) representing the output power of the filter 
and its white noise gain is optimized. 

13. The method of claim 12 wherein the relative weight of the 
output power and the "white noise gain" in the cost functional 
(J) is adjustable (using a parameter d) . 

14. The method of claim 12 wherein the cost functional (J) is 
minimized using the approximation that the sum, weighted by 
window functions, of the output of adjacent filters of the M 
filters is equal when applied to the same signal in time regions 
where said window functions overlap. 

15. The method of claim 1 comprising the further step of 
dividing filter coefficients into a fixed part and a adaptive 
part . 

16. The method of claim 1 comprising the further step of 
splitting the signals into at least two frequency bands, 
filtering at least one of said at least two frequency band, and, 
before further processing, recombining said frequency bands to 
regain noise attenuated seismic signals. 

17. The method of claim 1 wherein the discrete noisy signals 
used as input are recordings from individual seismic sensors 
prior to any group forming techniques. 



discrete seismic signals; and transmitting said signals as input 
to the filter. 

19. A method of processing seismic signals to generate a 
representation of a subterranean formation wherein noise 
present in said signals is attenuated by filtering said signals 
using the steps of - receiving signals using a plurality of 
receivers; 

- determining propagation characteristics of said signals with 
respect to receiver locations; and 

- filtering received signals using an at least partially 
adaptive filter such that signals having propagation 
characteristics other than said determined propagation 
characteristics are attenuated; 

said filtering step comprising the steps of 

- defining at least two independent sets of conditions with a 
first set defining a desired response and a second set defining 
said propagation characteristics of signals to be preserved ; 
and 

- adapting filter coefficients of said filter subject to said at 
least two independent sets of conditions so as to optimize the 
filter output for signals with propagation characteristics other 
than said determined propagation characteristics. 
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